#include <iostream>
#include <opencv2/core.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>

int main(int argc, char **argv) {
    cv::Mat mat = cv::imread("/home/xlll/Downloads/opencv/samples/data/lena.jpg", cv::IMREAD_GRAYSCALE);
    if (mat.empty()) {
        std::cout << "Failed to read image!" << std::endl;
        return EXIT_FAILURE;
    }
    
    /*
    cv::Mat kx = (cv::Mat_<float>(1, 2) << -1, 1);
    cv::Mat ky = (cv::Mat_<float>(2, 1) << -1, 1);
    cv::Mat gradx, grady;
    cv::filter2D(mat, gradx, CV_32FC1, kx, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::filter2D(mat, grady, CV_32FC1, ky, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    
    cv::Mat mags;
    cv::magnitude(gradx, grady, mags);
    cv::Mat angles = cv::Mat(mags.size(), CV_32FC1);
    for (int i = 0; i < gradx.rows; i++) {
        float* pAngle = angles.ptr<float>(i);
        float* pGradx = gradx.ptr<float>(i);
        float* pGrady = grady.ptr<float>(i);
        for (int j = 0; j < gradx.cols; j++) {
            pAngle[j] = std::atan2(pGrady[j], pGradx[j]) * 180. / CV_PI;
        }
    }
    cv::Mat gradx_8, grady_8, mags_8, angles_8;
    gradx.convertTo(gradx_8, CV_8UC1);
    grady.convertTo(grady_8, CV_8UC1);
    mags.convertTo(mags_8, CV_8UC1);
    angles.convertTo(angles_8, CV_8UC1);
    cv::imshow("mat", mat);
    cv::imshow("gradx", gradx_8);
    cv::imshow("grady", grady_8);
    cv::imshow("magnitude", mags_8);
    cv::imshow("angles", angles_8);
    */
    
    /*
    // Roberts
    cv::Mat k45 = (cv::Mat_<float>(2, 2) << 0, -1, 1, 0);
    cv::Mat k135 = (cv::Mat_<float>(2, 2) << -1, 0, 0, 1);
    cv::Mat grad45, grad135;
    cv::filter2D(mat, grad45, CV_8UC1, k45, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::filter2D(mat, grad135, CV_8UC1, k135, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::imshow("mat", mat);
    cv::imshow("grad45", grad45);
    cv::imshow("grad135", grad135);
    */
    
    /*
    // Prewitt
    cv::Mat kx = (cv::Mat_<float>(3, 3) << -1, 0, 1, -1, 0, 1, -1, 0, 1);
    cv::Mat ky = (cv::Mat_<float>(3, 3) << -1, -1, -1, 0, 0, 0, 1, 1, 1);
    cv::Mat k45 = (cv::Mat_<float>(3, 3) << 0, -1, -1, 1, 0, -1, 1, 1, 0);
    cv::Mat k135 = (cv::Mat_<float>(3, 3) << -1, -1, 0, -1, 0, 1, 0, 1, 1);
    cv::Mat gradx, grady, grad45, grad135;
    cv::filter2D(mat, gradx, CV_32FC1, kx, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::filter2D(mat, grady, CV_32FC1, ky, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::filter2D(mat, grad45, CV_32FC1, k45, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::filter2D(mat, grad135, CV_32FC1, k135, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::Mat mags;
    cv::magnitude(gradx, grady, mags);
    cv::Mat angles = cv::Mat(mags.size(), CV_32FC1);
    for (int i = 0; i < gradx.rows; i++) {
        float* pAngle = angles.ptr<float>(i);
        float* pGradx = gradx.ptr<float>(i);
        float* pGrady = grady.ptr<float>(i);
        for (int j = 0; j < gradx.cols; j++) {
            pAngle[j] = std::atan2(pGrady[j], pGradx[j]) * 180 / CV_PI;
        }
    }
    cv::Mat gradx_8, grady_8, mags_8, angles_8, grad45_8, grad135_8;
    gradx.convertTo(gradx_8, CV_8UC1);
    grady.convertTo(grady_8, CV_8UC1);
    mags.convertTo(mags_8, CV_8UC1);
    angles.convertTo(angles_8, CV_8UC1);
    grad45.convertTo(grad45_8, CV_8UC1);
    grad135.convertTo(grad135_8, CV_8UC1);
    cv::imshow("mat", mat);
    cv::imshow("gradx", gradx_8);
    cv::imshow("grady", grady_8);
    cv::imshow("magnitude", mags_8);
    cv::imshow("angles", angles_8);
    cv::imshow("grad45", grad45_8);
    cv::imshow("grad135", grad135_8);
    */
    
    /*
    // Sobel
    cv::Mat kx_x, kx_y, kx;
    cv::getDerivKernels(kx_x, kx_y, 1, 0, 3, false, CV_32FC1);
    kx = kx_y * kx_x.t();
    std::cout << "kx_x = \n" << kx_x << std::endl;
    std::cout << "kx_y = \n" << kx_y << std::endl;
    std::cout << "kx = \n" << kx << std::endl;
    cv::Mat gradx1, gradx2, gradx3;
    cv::sepFilter2D(mat, gradx1, CV_32FC1, kx_x, kx_y, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::filter2D(mat, gradx2, CV_32FC1, kx, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::Sobel(mat, gradx3, CV_32FC1, 1, 0, 3, 1.0, 0.0, cv::BORDER_REFLECT101);
    cv::Mat diffx1, diffx2;
    diffx1 = gradx1 != gradx2;
    diffx2 = gradx1 != gradx3;
    std::vector<cv::Point> nonzerox1, nonzerox2;
    cv::findNonZero(diffx1, nonzerox1);
    cv::findNonZero(diffx2, nonzerox2);
    std::cout << "nonzerox1 size = " << nonzerox1.size() << std::endl;
    std::cout << "nonzerox2 size = " << nonzerox2.size() << std::endl;
    
    cv::Mat ky_x, ky_y, ky;
    cv::getDerivKernels(ky_x, ky_y, 0, 1, 3, false, CV_32FC1);
    ky = ky_y * ky_x.t();
    std::cout << "ky_x = \n" << ky_x << std::endl;
    std::cout << "ky_y = \n" << ky_y << std::endl;
    std::cout << "ky = \n" << ky << std::endl;
    cv::Mat grady1, grady2, grady3;
    cv::sepFilter2D(mat, grady1, CV_32FC1, ky_x, ky_y, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::filter2D(mat, grady2, CV_32FC1, ky, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::Sobel(mat, grady3, CV_32FC1, 0, 1, 3, 1.0, 0.0, cv::BORDER_REFLECT101);
    cv::Mat diffy1, diffy2;
    diffy1 = grady1 != grady2;
    diffy2 = grady1 != grady3;
    std::vector<cv::Point> nonzeroy1, nonzeroy2;
    cv::findNonZero(diffy1, nonzeroy1);
    cv::findNonZero(diffy2, nonzeroy2);
    std::cout << "nonzeroy1 size = " << nonzeroy1.size() << std::endl;
    std::cout << "nonzeroy2 size = " << nonzeroy2.size() << std::endl;
    
    cv::Mat mags;
    cv::magnitude(gradx1, grady1, mags);
    cv::Mat angles = cv::Mat(mags.size(), CV_32FC1);
    for (int i = 0; i < gradx1.rows; i++) {
        float* pAngle = angles.ptr<float>(i);
        float* pGradx = gradx1.ptr<float>(i);
        float* pGrady = grady1.ptr<float>(i);
        for (int j = 0; j < gradx1.cols; j++) {
            pAngle[j] = std::atan2(pGrady[j], pGradx[j]) * 180 / CV_PI;
        }
    }
    cv::Mat gradx_8, grady_8, mags_8, angles_8;
    gradx1.convertTo(gradx_8, CV_8UC1);
    grady1.convertTo(grady_8, CV_8UC1);
    mags.convertTo(mags_8, CV_8UC1);
    angles.convertTo(angles_8, CV_8UC1);
    cv::imshow("mat", mat);
    cv::imshow("gradx", gradx_8);
    cv::imshow("grady", grady_8);
    cv::imshow("magnitude", mags_8);
    cv::imshow("angles", angles_8);
    
    cv::Mat k45 = (cv::Mat_<float>(3, 3) << 0, -1, -2, 1, 0, -1, 2, 1, 0);
    cv::Mat k135;
    cv::flip(k45, k135, 1);
    std::cout << "k45 = \n" << k45 << std::endl;
    std::cout << "k135 = \n" << k135 << std::endl;
    cv::Mat grad45, grad135;
    cv::filter2D(mat, grad45, CV_32FC1, k45, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::filter2D(mat, grad135, CV_32FC1, k135, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::Mat grad45_8, grad135_8;
    grad45.convertTo(grad45_8, CV_8UC1);
    grad135.convertTo(grad135_8, CV_8UC1);
    cv::imshow("grad45", grad45_8);
    cv::imshow("grad135", grad135_8);
    */
    
    //Scharr
    cv::Mat kx_x, kx_y, kx;
    cv::getDerivKernels(kx_x, kx_y, 1, 0, cv::FILTER_SCHARR, false, CV_32FC1);
    kx = kx_y * kx_x.t();
    std::cout << "kx_x = \n" << kx_x << std::endl;
    std::cout << "kx_y = \n" << kx_y << std::endl;
    std::cout << "kx = \n" << kx << std::endl;
    cv::Mat gradx1, gradx2, gradx3;
    cv::sepFilter2D(mat, gradx1, CV_32FC1, kx_x, kx_y, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::filter2D(mat, gradx2, CV_32FC1, kx, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::Sobel(mat, gradx3, CV_32FC1, 1, 0, cv::FILTER_SCHARR, 1.0, 0.0, cv::BORDER_REFLECT101);
    cv::Mat diffx1, diffx2;
    diffx1 = gradx1 != gradx2;
    diffx2 = gradx1 != gradx3;
    std::vector<cv::Point> nonzerox1, nonzerox2;
    cv::findNonZero(diffx1, nonzerox1);
    cv::findNonZero(diffx2, nonzerox2);
    std::cout << "nonzerox1 size = " << nonzerox1.size() << std::endl;
    std::cout << "nonzerox2 size = " << nonzerox2.size() << std::endl;
    
    cv::Mat ky_x, ky_y, ky;
    cv::getDerivKernels(ky_x, ky_y, 0, 1, cv::FILTER_SCHARR, false, CV_32FC1);
    ky = ky_y * ky_x.t();
    std::cout << "ky_x = \n" << ky_x << std::endl;
    std::cout << "ky_y = \n" << ky_y << std::endl;
    std::cout << "ky = \n" << ky << std::endl;
    cv::Mat grady1, grady2, grady3;
    cv::sepFilter2D(mat, grady1, CV_32FC1, ky_x, ky_y, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::filter2D(mat, grady2, CV_32FC1, ky, cv::Point(-1, -1), 0.0, cv::BORDER_REFLECT101);
    cv::Sobel(mat, grady3, CV_32FC1, 0, 1, cv::FILTER_SCHARR, 1.0, 0.0, cv::BORDER_REFLECT101);
    cv::Mat diffy1, diffy2;
    diffy1 = grady1 != grady2;
    diffy2 = grady1 != grady3;
    std::vector<cv::Point> nonzeroy1, nonzeroy2;
    cv::findNonZero(diffy1, nonzeroy1);
    cv::findNonZero(diffy2, nonzeroy2);
    std::cout << "nonzeroy1 size = " << nonzeroy1.size() << std::endl;
    std::cout << "nonzeroy2 size = " << nonzeroy2.size() << std::endl;
    
    cv::Mat mags;
    cv::magnitude(gradx1, grady1, mags);
    cv::Mat angles = cv::Mat(mags.size(), CV_32FC1);
    for (int i = 0; i < gradx1.rows; i++) {
        float* pAngle = angles.ptr<float>(i);
        float* pGradx = gradx1.ptr<float>(i);
        float* pGrady = grady1.ptr<float>(i);
        for (int j = 0; j < gradx1.cols; j++) {
            pAngle[j] = std::atan2(pGrady[j], pGradx[j]) * 180 / CV_PI;
        }
    }
    cv::Mat gradx_8, grady_8, mags_8, angles_8;
    gradx1.convertTo(gradx_8, CV_8UC1);
    grady1.convertTo(grady_8, CV_8UC1);
    mags.convertTo(mags_8, CV_8UC1);
    angles.convertTo(angles_8, CV_8UC1);
    cv::imshow("mat", mat);
    cv::imshow("gradx", gradx_8);
    cv::imshow("grady", grady_8);
    cv::imshow("magnitude", mags_8);
    cv::imshow("angles", angles_8);
    
    cv::waitKey(0);
    return EXIT_SUCCESS;
}
